#ifndef ERK_HPP
#define ERK_HPP

#include <functional>
#include <iostream>

#include "RK.hpp"

class ERK : public RK {
public:
  ERK(int dim, double dt, double t0, double t1) : RK(dim, dt, t0, t1) {}

  void solve(const Eigen::VectorXd &u0,
             const std::function<Eigen::VectorXd(Eigen::VectorXd, double)> &f) {
    if (u0.size() != m_dim) {
      std::cerr << "invalid input size." << std::endl;
    }
    m_u.row(0) = u0;
    Eigen::MatrixXd Y(m_s, m_dim);
    for (int i = 1; i < m_u.rows(); ++i) {
      
      double t = m_t0 + (i - 1) * m_dt;
      Y.setZero();
      for (int j = 0; j < m_s; ++j) {
        Eigen::VectorXd tmp(m_u.row(i - 1));
        for (int k = 0; k < j; ++k) {
          tmp += m_dt * m_A(j, k) * Y.row(k);
        }
        Y.row(j) = f(tmp, t + m_c(j) * m_dt);
      }
      m_u.row(i) = m_u.row(i - 1);
      for (int j = 0; j < m_s; ++j) {
        m_u.row(i) += m_dt * m_b(j) * Y.row(j);
      }
    }
  }
};

#endif